/**
  CNOK project, Anyang Normal University, IMP-CAS
  \class TABound
  \brief To calculate bound state radial wavefunction of a valence nucleon in
  central potential (Vc+VN+VLS) of a nucleus.
  \author SUN Yazhou, asia.rabbit@163.com
  \since 2020/07/09
  \date Last modified: 2022/02/07 by SUN Yazhou
  \copyright 2020-2023 SUN Yazhou
  \copyright CNOK project, Anyang Normal University, IMP-CAS
*/

#include "config.h"
#ifdef CATCH_CNOK
#include "catch.h"
#endif
#include <algorithm>
#include <iostream>
#include <cstring>
#include <fstream>
#include "TABound.h"
#include "TAMath.h"
#include "TAMinimize.h"
#include "TAInterpolate.h"
#include "TAEqSetSolver.h"
#include "TAIntegrate.h"
#include "TAAME.h"
#include "TAParaManager.h"
#include "TAEqSetSolver.h"

using std::max;
using std::ofstream;
using std::cout;
using std::endl;
using std::setprecision;
using std::fixed;

#define polint TAInterpolateD::PolyInter
#define romb TAIntegrate<>::Romberg
#define simp TAIntegrate<>::Simpson
#define sqr TAMath::sqr

// physical constants for calculation of the potential
static const double HBARC = TAMath::hbarc();
static const double RP2 = sqr(HBARC/TAMath::PionPMMass()); // pion compton wavelength^2
static const double e2 = HBARC*TAMath::FineStructureConstant(); // e^2/(4pi*epsilon0)
static const TAAME *AME = TAAME::Instance();
// static const TAParaManager *PARA = TAParaManager::Instance();
static const char LL[] = {'s', 'p', 'd', 'f', 'g', 'h', 'i', 'j', 'k', 'l'};
static int NBOUND = 0; // bookkeeping calling of TABound::Bound()

double TABound::ZERO = 1.e-6; // in fm, the x1, physical origin to nuclear systems
double TABound::INFTY = 50.; // in fm, the x2, physical infinity to nuclear systems

TABound::TABound(const TAYaml &v) : TAODESolver(), fSolved(false), fxx(nullptr),
fyy2(nullptr), fcnt(){
  // set integration range //
  if(v["ZEROR"]) ZERO = v["ZEROR"].as<double>();
  if(v["INFTYR"]) INFTY = v["INFTYR"].as<double>();
  // (c+v) is P (projectile) //
  fZc = v["Zc"].as<int>();
  fAc = v["Ac"].as<int>();
  fZv = v["ZP"].as<int>() - fZc;
  fAv = v["AP"].as<int>() - fAc;
  // the quantum state of the valence nucleon //
  n = v["n"].as<int>();
  l = v["l"].as<int>();
  j = v["j"].as<double>();
  // assign the potential parameters //
  // a scaling factor 2\mu/hbar^2 //
  const double acM = AME->massNUC(fZc, fAc);
  const double avM = AME->massNUC(fZv, fAv);
  const double mu = acM*avM/(acM+avM); // the reduced mass in MeV
  const double m2 = f2MuOverHbar2 = 2.*mu/(HBARC*HBARC);

  // configure the potential //
  fEref = 0.; fRref = 0.;
  if(v["Eref"]) fEref = -v["Eref"].as<double>()*f2MuOverHbar2;
  if(v["RHFrms"]) fRref = v["RHFrms"].as<double>()*sqrt((fAc+fAv)/fAc);
  const double DR0 = (2*n+l-j*0.25)*0.5; // emperical formula
  const double R0 = fRref ? fRref*1.2 + DR0 : 1.6*pow(fAc, 1./3) + DR0;
  // cout << "R0: " << R0 << endl; getchar();
  fV0 = (v["V0"] ? v["V0"].as<double>() : -50.) * m2;
  fR0 = v["r0"] ? v["r0"].as<double>() : R0;
  fRS = v["rS"] ? v["rS"].as<double>() : R0;
  fRC = v["rC"] ? v["rC"].as<double>() : R0;
  fA0 = v["a"].as<double>();
  fVS = v["VS"].as<double>() * m2;
  fAS = v["aS"].as<double>();
  fFindOpt = 0;
  if(v["findOpt"]) fFindOpt = v["findOpt"].as<int>();
  // the orbit info //
  char orbit[256];
  sprintf(orbit, "%d%c%d/2_%s_%s", n, LL[l], int(j*2),
    AME->name(fZc, fAc), AME->name(fZv, fAv));
  fOrbit = orbit;

  fCoel = l*(l+1.);
  fCoeVS = fVS*RP2/fAS * (j*(j+1.)-fCoel-0.75)/2.;
  fCoeCoul = fZc*fZv*e2 * m2;

  BoundIter(); // solve the bound state
} // end the constructor

TABound::~TABound(){
  if(fxx) delete [] fxx;
  if(fyy2) delete [] fyy2;
} // end the dtor

// self-adaptive BoundIter according to availabilities of fRref and fRref
void TABound::BoundIter(){
  if(0. == fEref && 0. == fRref) return Bound();
  if(0. != fEref && 0. != fRref){
    if(0 == fFindOpt) return BoundIterNewton();
    else return BoundIterER();
  } // end if
  if(0. != fEref) return BoundIterE();
  TAException::Error("TABound", "BoundIter: fRref is provided while fEref is not.\
\nIn which case we do not serve to only reproduce Rrms.");
} // end member function BoundIter

void TABound::BoundIterE(){
  if(0. == fEref) TAException::Error("TABound", "BoundIterE: fEref is 0.!");
  // optimize for E //
  static const int NV = 1;
  auto f = [=](int n, const double *x, double *fvec){
    fV0 = x[0]; fSolved = false; Bound(); // reinitiate Bound()
    fvec[0] = fabs(fE - fEref)/fEref; // + fabs(fRrms - fRref)/fRref*0.3;
  };
  double tolx, tolf;
  TAEqSetSolver::GetTOL(tolx, tolf);
  TAEqSetSolver::SetTOL(1.e-2, 1.e-2);
  double x[NV] = {fV0};
  int check = TAEqSetSolver::Newton(NV, x, f);
  TAEqSetSolver::SetTOL(tolx, tolf); // recover tol
  if(check) TAException::Warn("TABound", "BoundIterE: Check for mere trapped in a minimum instead of hitting 0.");
  cout << "______####_____fE:    " << fE/f2MuOverHbar2 << " MeV" << endl;
  cout << "______####_____fEref: " << fEref/f2MuOverHbar2 << " MeV" << endl;
  cout << "\033[33;1m_____________fV0 = " << fV0 / f2MuOverHbar2 << " MeV" << endl;
  cout << "\033[33;1m_____________fR0 = " << fR0 << " fm\033[0m" << endl;
  cout << endl << "V0: " << fV0 / f2MuOverHbar2 << endl;
  cout << "Totally " << NBOUND << " call(s) of Bound()." << endl;
  cout << "r0: " << fR0 << endl;
  cout << "rS: " << fR0 << endl;
  cout << "rC: " << fR0 << endl;
  // getchar(); // exit(1);
} // end member function BoundIterE
void TABound::BoundIterR(){
  if(0. == fRref) TAException::Error("TABound", "BoundIterR: fRref is 0.!");
  // optimize for E //
  static const int NV = 1;
  auto f = [=](int n, const double *x, double *fvec){
    fR0 = fRS = fRC = x[0]; fSolved = false; Bound(); // reinitiate Bound()
    fvec[0] = fabs(fRrms - fRref)/fRref;
  };
  double tolx, tolf;
  TAEqSetSolver::GetTOL(tolx, tolf);
  TAEqSetSolver::SetTOL(1.e-3, 1.e-3);
  double x[NV] = {fR0};
  int check = TAEqSetSolver::Newton(NV, x, f);
  TAEqSetSolver::SetTOL(tolx, tolf); // recover tol
  if(check) TAException::Warn("TABound", "BoundIterE: Check for mere trapped in a minimum instead of hitting 0.");
  cout << "______####_____fRrms: " << fRrms << " fm" << endl;
  cout << "______####_____fRref: " << fRref << " fm" << endl;
  cout << "\033[33;1mfV0 = " << fV0 / f2MuOverHbar2 << " MeV" << endl;
  cout << "\033[33;1mfR0 = " << fR0 << " fm\033[0m" << endl;
  cout << "Totally " << NBOUND << " call(s) of Bound()." << endl;
  cout << endl << "V0: " << fV0 / f2MuOverHbar2 << endl;
  cout << "r0: " << fR0 << endl;
  cout << "rS: " << fR0 << endl;
  cout << "rC: " << fR0 << endl;
  // getchar(); // exit(1);
} // end member function BoundIterR

void TABound::BoundIterER(){
  if(0. == fEref || 0. == fRref) TAException::Error("TABound", "BoundIterER: fEref or fEref is 0.!");
  cout << "______####_____fEref: " << fEref/f2MuOverHbar2 << " MeV" << endl;
  cout << "______####_____fRref: " << fRref << " fm" << endl;
  const double R0 = 1.25*pow(fAc, 1./3.), DR0 = R0*0.25;
  double r0m = 0.;
  //----------------- PLAIN SEARCH ----------------------//
  // double drmin = 1.e200;
  // const int nn = 20;
  // for(int i = 0; i <= nn; i++){
  //   const double r0 = R0+DR0*(2.*i/double(nn)-1.);
  //   fR0 = fRS = fRC = r0;
  //   BoundIterE();
  //   if(fabs(fRrms-fRref) < drmin){
  //     drmin = fabs(fRrms-fRref);
  //     r0m = r0;
  //   } // end if
  // } // end for over i
  //----------------- EXQUISITE SEARCH -------------------//
  auto f = [=](double r0){ fR0 = fRS = fRC = r0; BoundIterE(); return fabs(fRrms-fRref); };
  // TAMinimize<>::Golden(R0-DR0, R0, R0+DR0, f, r0m, 1.e-2);
  TAMinimize<>::Brent(R0-DR0, R0, R0+DR0, f, r0m, 1.e-2);
  //------------------ DONE DONE DONE ------------------ -//

  fR0 = fRS = fRC = r0m; BoundIterE();
  // BoundIterR();
  cout << "______####_____fE:    " << fE/f2MuOverHbar2 << " MeV" << endl;
  cout << "______####_____fEref: " << fEref/f2MuOverHbar2 << " MeV" << endl;
  cout << "______####_____fRrms: " << fRrms << " fm" << endl;
  cout << "______####_____fRref: " << fRref << " fm" << endl;
  cout << "R0: " << R0 << " R0MIN: " << R0 - DR0 << " R0MAX: " << R0 + DR0 << endl;
  cout << "\033[33;1mfV0 = " << fV0 / f2MuOverHbar2 << " MeV" << endl;
  cout << "\033[33;1mfR0 = " << fR0 << " fm\033[0m" << endl;
  cout << "Totally " << NBOUND << " call(s) of Bound()." << endl;
  cout << endl << "V0: " << fV0 / f2MuOverHbar2 << endl;
  cout << "r0: " << fR0 << endl;
  cout << "rS: " << fR0 << endl;
  cout << "rC: " << fR0 << endl;
  getchar();
  // getchar(); exit(1);
} // end member function BoundIterER
void TABound::BoundIterNewton(){
  if(0. == fEref || 0. == fRref) TAException::Error("TABound", "BoundIterER1: fEref or fEref is 0.!");
  // optimize for E and Rrms //
  static const int NV = 2;
  auto f = [=](int n, const double *x, double *fvec){
    fV0 = x[0]; fR0 = fRS = fRC = x[1];
    fSolved = false; Bound(); // reinitiate Bound()
    fvec[0] = fabs(fE - fEref)/fEref;
    fvec[1] = fabs(fRrms - fRref)/fRref;
  };
  double tolx, tolf;
  TAEqSetSolver::GetTOL(tolx, tolf);
  TAEqSetSolver::SetTOL(1.e-3, 1.e-3);
  double x[NV] = {fV0, fR0};
  int check = TAEqSetSolver::Newton(NV, x, f);
  TAEqSetSolver::SetTOL(tolx, tolf); // recover tol
  if(check) TAException::Warn("TABound", "BoundIterER1: Check for mere trapped in a minimum instead of hitting 0.");
  cout << "______####_____fE:    " << fE/f2MuOverHbar2 << " MeV" << endl;
  cout << "______####_____fEref: " << fEref/f2MuOverHbar2 << " MeV" << endl;
  cout << "______####_____fRrms: " << fRrms << " fm" << endl;
  cout << "______####_____fRref: " << fRref << " fm" << endl;
  cout << "\033[33;1mfV0 = " << fV0 / f2MuOverHbar2 << " MeV" << endl;
  cout << "\033[33;1mfR0 = " << fR0 << " fm\033[0m" << endl;
  cout << "Totally " << NBOUND << " call(s) of Bound()." << endl;
  cout << endl << "V0: " << fV0 / f2MuOverHbar2 << endl;
  cout << "r0: " << fR0 << endl;
  cout << "rS: " << fR0 << endl;
  cout << "rC: " << fR0 << endl;
  // getchar();
} // end member function BoundIterNewton

void TABound::derivs(double r, const double *y, double *dydx){
  dydx[0] = y[1];
  dydx[1] = (VT(r)-fE)*y[0];
} // end member function derivs

///< \retval V(r) = V0+VSO+VC+VL (VL is the centrifugal barrier)
double TABound::V(double r) const{
  const double v0 = fV0/(1.+exp((r-fR0)/fA0)); // the central NN potential
  // the spin-orbin potential //
  double vs = exp((r-fRS)/fAS); // vs = -vs0*r^{pi}_{compton}^2/r*dws/dr
  vs = fCoeVS * vs/(sqr(1.+vs)*r);
  // the Coulomb potential //
  double vc = fCoeCoul;
  if(r >= fRC) vc /= r; else vc *= (3.-sqr(r/fRC)) / (2.*fRC);

  return v0 + vs + vc;
} // end member function V
double TABound::VT(double r) const{
  double v = V(r);
  if(l){
    if(0. == r) TAException::Error("TABound", "VT: r is 0 for nonzero l");
    v += fCoel/(r*r);
  } // end if(l)
  return v;
} // end member function VT

// #define DEBUG

void TABound::Bound(){
  if(fSolved) return;
  fSolved = true; // has been called
  NBOUND++;

  cout << endl << "Solving bound state (nlj) " << n << LL[l] << int(j*2.) << "/2 (c+v) ";
  cout << AME->name(fZc, fAc) << "+" << AME->name(fZv, fAv) << "..." << endl;

  SetEPS(1.e-6); // relative error at each step of RKstepper
  SetInitialStep(1.e-5);
  const double r0 = ZERO, rInf = INFTY;
  SetMinStep(1.e-30);
  // (int maxStepCount, double mindxsave, int nvar)
  // the mindxsave is so that the x-range covers the whole interval //
  const int NMAX = 5000; // maximum steps allowed in saving the ODE result
  SetSave(NMAX, (rInf-r0)/NMAX, 2); // (nmax, h1, nvar).
  // n = 2; l = 2;
  // the initial guess for the eigenenergy
  // static const double EREF = -4.99*f2MuOverHbar2;
  // calculate min of V for later use
  // draw VT(r) //
  ofstream foutv("vis/VT.txt"); // the total potential VT(r)
  const int n_ = 500; const double dr_ = (10.)/n_;
  for(int i = n_; i--;){
    const double r = dr_*i+0.6;
    foutv << r << " " << VT(r)/f2MuOverHbar2 << endl;
  } // end for over i
  foutv.close();
  double rmin = 0., vmin = 0.; // VT(rmin)=vmin
  const auto& f = [=](double r){ return VT(r); };
  vmin = TAMinimize<>::Brent(r0, 2., 7., f, rmin);
  // V+T=E, T=E-V>0, so V<E != empty_set is mandatory for physical solutions //
  if(vmin >= fEref) TAExcp::Error("TABound",
    "Bound: fEref(%f) > vmin(%f) is not physical", fEref, vmin);
  if(0. == fEref) fE = 0.5*vmin; // not assigned to Eref, then assign it to 0.5vmin
  else fE = fEref;
  // cout << "VT_min(" << rmin << ")/2 = " << fE << endl; //- DEBUG
  cout << "Firstly, fE = " << fE/f2MuOverHbar2 << " MeV" << endl;
  double rm = 0.; // the turning point of total effective potential
  double dE = 1.e200;
  fxx = new double[NMAX*2]{}; // the concatenated solution: x
  fyy = new double[NMAX*2]{}; // the concatenated solution: y
  fyy2 = new double[NMAX*2]{}; // the concatenated solution: y^2
  // arrays to store the solutions //
  double ri[NMAX]{}, ui[NMAX]{}; // the interior solution: r0->rm
  double re[NMAX]{}, ue[NMAX]{}; // the exterior solution: rInf->rm
  double r2y2[NMAX*2]{}; // calculate r_rms, r2y2: (r*u(r))^2
  // iteration to find the optimal eigenenergy //
  cout << setprecision(6) << fixed; // DEBUG
  char name1[] = "vis/rRl.txt";
  char name2[] = "vis/rRl2.txt";
  char name3[] = "vis/Rl.txt";
  int niter = 0; // too many iterations indicates accumulation of roundoff, and catastrophic divergence
  while(fabs(dE) > GetEPS()){ // usually E has the same accuracy with EPS
    if(niter++ > 40.) TAException::Warn("TABound", "Bound: too many iterations occurred!");
    ofstream fout1(name1), fout2(name2), fout3(name3);
    // kinetic energy is postive at the inner part of the potential //
    if(fE < V(r0)) fE = V(r0) + fabs(V(r0)*0.1); // an emperical try
    // calculate the turning point rm as the first point from infinity where f(rm) = 0 //
    const double dr = rInf/5000;
    int i = NMAX;
    for(; i > 0; i--) if(VT(dr*i)-fE < 0.){ rm = dr*i; break; }
    if(i == 0) rm = 1.; // expediently set to 1
    // cout << "rm: " << rm << "fm" << endl; getchar(); // DEBUG
#ifdef DEBUG
    cout << "VT(" << rm << ")-" << fE << ": " << VT(rm)-fE << endl; //- DEBUG
#endif
    //-- solve the ODE --//
    // integrate from r0 to rm //
    const double k = sqrt(fE-V(r0));
    const double x = r0*k, x2 = x*x, x4 = x2*x2, x6 = x4*x2;
    double l3 = (2.*l+3.), l5 = l3*(2.*l+5.), l7 = l5*(2.*l+7.);
    l3 = x2/(2.*l3); l5 = x4/(8.*l5); l7 = x6/(48.*l7);
    fy[0] = x*(1.-l3+l5-l7); // accuracy: O(x^8)*x
    fy[1] = k*(l+1.)-l3*(l+3.)+l5*(l+5.)-l7*(l+7.); // accuracy: O(x^8)
    ODEIntegrator(r0, rm);
    const int cnti = GetSolution(ri, ui); // extract the solution
    const double li = fy[1]/fy[0];
    // integrate from rInf to rm //
    const double beta = sqrt(fabs(fE));
    fy[0] = 1.e-6; // exp(-beta*rInf)
    fy[1] = -beta*fy[0]; // -beta*
    ODEIntegrator(rInf, rm);
    const int cnte = GetSolution(re, ue);
    const double le = fy[1]/fy[0];
    // concatenate the solution //
    for(int i = 0; i < cnti; i++){
      fxx[i] = ri[i];
      fyy[i] = ui[i];
      fyy2[i] = sqr(fyy[i]);
    } // end for
    const double sca = ui[cnti-1]/ue[cnte-1];
    for(int i = 1; i < cnte; i++){
      const int i0 = cnti+i-1, i1 = cnte-i-1; // r=rm is dropped
      fxx[i0] = re[i1];
      fyy[i0] = ue[i1] * sca;
      fyy2[i0] = sqr(fyy[i0]);
    } // end for
    const int cnt = cnti+cnte-1; fcnt = cnt;
    const int nodes = TAMath::countNodes(fyy, cnt);
    // normalize the function, and calculate r_{rms}: sqrt(<psi|r^2|psi>) //
    // XXXXXXXXXXX NEVER USE SIMPSON'S RULE HERE, FOR fx IS _NOT_ EQUALLY SPACED XXXXXXXXXXXXX
    const double cc = pow(romb(cnt, fxx, fyy2, 1.e-3), -0.5); // +fxx[0]*fyy2[0]/2.
    for(int i = cnt; i--;){
      fyy[i] *= cc;
      fyy2[i] = sqr(fyy[i]);
      r2y2[i] = fyy2[i]*sqr(fxx[i]);
    } // end for over i
    // XXXXXXXXXXX NEVER USE SIMPSON'S RULE HERE, FOR fx IS _NOT_ EQUALLY SPACED XXXXXXXXXXXXX
    fRrms = sqrt(romb(cnt, fxx, r2y2, 1.e-3)); // +fxx[0]*r2y2[0]/2.
    cout << "_______rms: " << fRrms << " fm" << endl;
    //////////// DEBUG DEBUG DEBUG DEBUG DEBUG DEBUG DEBUG DEBUG DEBUG //
#ifdef DEBUG
    cout << "li: " << li << " \tle: " << le << endl; // DEBUG
    cout << "li-le: " << li-le << endl; // DEBUG
#endif
    for(int i = 0; i < cnt; i++) fout1 << fxx[i] << " " << fyy[i] << endl; // DEBUG
    for(int i = 0; i < cnt; i++) fout2 << fxx[i] << " " << fyy2[i] << endl; // DEBUG
    for(int i = 0; i < cnt; i++) fout3 << fxx[i] << " " << fyy[i]/fxx[i] << endl; // DEBUG
#ifdef DEBUG
    cout << "Rl written to file" << endl; // DEBUG
    cout << "nodes: " << nodes << endl; // DEBUG
    //////////// DEBUG DEBUG DEBUG DEBUG DEBUG DEBUG DEBUG DEBUG DEBUG //
    cout << "______fE_before: " << fE/f2MuOverHbar2 << " MeV" << endl; // DEBUG
#endif
    if(nodes > n) fE *= 1.1;
    if(nodes < n) fE *= 0.9;
    if(nodes == n){ // calculate and apply the energy correction //
      dE = (li-le)*fyy2[cnti-1];
#ifdef DEBUG
      cout << "initial dE: " << dE << endl; // DEBUG
#endif
      // fE < 0 to ensure bound states
      // fE > vmin since E=V+T, E_min=V_min, as T_min=0
      if(fabs(dE/fE) > 1.e-20) while(fE+dE > 0. || fE+dE < vmin) dE *= 0.5;
      fE += dE;
    } // end else
#ifdef DEBUG
    // cout << "Eref: " << EREF << endl; // DEBUG
    // cout << "______fE_after: " << fE << endl; // DEBUG
#endif
    cout << "______fE_after: " << fE/f2MuOverHbar2 << " MeV" << endl; // DEBUG
    cout << "dE: " << dE/f2MuOverHbar2 << " MeV" << endl; // DEBUG
    // getchar(); // DEBUG
    fout1.close(); fout2.close(); fout3.close();
  } // end while
  cout << "\033[32;1m______E: " << fE/f2MuOverHbar2 << " MeV\033[0m" << endl; // DEBUG
  cout << "\033[34;1m_______rms: " << fRrms << " fm\033[0m" << endl; // DEBUG
  cout << endl << "...DONE. r*Rl(r) stored in \033[32;1m" << name1 << "\033[0m." << endl;
  cout << "|r*Rl(r)|^2 stored in \033[32;1m" << name2 << "\033[0m." << endl;
  cout << "Rl(r) stored in \033[32;1m" << name3 << "\033[0m." << endl << endl;
  cout << "Rl(r), rRl(r) and |rRl(r)|^2 can be inspected via `cd vis/`";
  cout << " and `root -l Rl.C`" << endl;
  cout << "(provided ROOT is installed)." << endl;
  cout << "\033[32;43;1m___BINGO___BINGO___BINGO\033[0m" << endl; // DEBUG
} // end member function Bound

/// \retval Note that Rl=u/r; u = Rl*r. This method returns Rl, i.e. u/r, NOT u
/// implement an interpolation
double TABound::u(double r){ // u=r*Rl
  Bound();

  static const int NP = 5; // polynomial interpolation: Nof points
  if(0. == r) r = 1.e-5;
  if(r < 0.) TAException::Error("TABound", "Rl: input r(%f) is minus.", r);
  if(r > fxx[fcnt-1]) TAException::Error("TABound",
    "Rl: input r(%f) out of bounds(r_max=%f).", r, fxx[fcnt-1]);
  const double rl = polint(fxx, fyy, fcnt, NP, r, nullptr, false);
  return fabs(rl) < 1.e-5 ? 0. : rl;
} // end member function u

double TABound::Rl(double r){
  if(0. == r) r = 1.e-5;
  return u(r)/r;
} // end member function Rl


#ifdef CATCH_CNOK
// test the minimizer Brent //
TEST_CASE("Brent minimizer", "[brent]"){
  double xmin = 0., ymin = 0.;
  double a = -1., b = 4.;
  const auto& f = [](double x){ return (x-1.)*(x+2.)*(x-3.); };
  ymin = TAMinimize<>::Brent(a, a, b, f, xmin, 1.e-8);
  double ex = (2.+sqrt(19.))/3., ey = (56.-sqrt(27436.))/27.;
  // cout << "xmin: " << xmin << " ymin: " << ymin << endl;
  CHECK(xmin == Approx(ex).epsilon(1.e-4));
  CHECK(ymin == Approx(ey).epsilon(1.e-4));
} // end TEST_CASE brent
#endif
